arXiv: 1507.05993v 1 [nlin.CD] 9 Jul 2015 


Local observability of state variables and parameters in nonlinear modeling 
quantified by delay reconstruction 

Ulrich Parlitz, 1,2, [2] Jan Schumann-Bischoff/^and Stefan Luther 1 ' 

^ Max Planck Institute for Dynamics and Self-Organization, Am Fafiberg 17, 37077 Gottingen, 

Germany 

2 )Institute for Nonlinear Dynamics, G eorg-August-Universitat Gottingen, Am Fafiberg 17, 37077 Gottingen, 
Germany 

(Dated: 23 July 2015) 

Features of the Jacobian matrix of the delay coordinates map are exploited for quantifying the robustness 
and reliability of state and parameter estimations for a given dynamical model using an observed time series. 
Relevant concepts of this approach are introduced and illustrated for discrete and continuous time systems 
employing a filtered Henon map and a Rossler system. 
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For many physical processes dynamical models 
(differential equations or iterated maps) are avail¬ 
able but often not all of their variables and pa¬ 
rameters are known or can be (easily) measured. 
In meteorology, for example, sophisticated large 
scale models exist, which have to be continuously 
adapted to the true temporal changes of temper¬ 
atures, wind speed, humidity, and other relevant 
physical quantities. To obtain a model that “fol¬ 
lows” reality, measured data have to be repeat¬ 
edly incorporated into the model. In geosciences 
this procedure is called data assimilation , but the 
task to track state variables and system parame¬ 
ters by means of estimation methods occurs also 
in other fields of physics and applications. How¬ 
ever, not all observables provide the information 
required to estimate a particular unknown quan¬ 
tity. In this article, we consider this problem of 
observability in the context of chaotic dynamics 
where sensitive dependance on initial conditions 
complicates any estimation method. A quanti¬ 
tative characterization of local observability em¬ 
ploying delay coordinates is used to answer the 
question where in state and parameter space esti¬ 
mation of a particular state variable or parameter 
is feasible and where not. 


I. INTRODUCTION 

To describe and forecast dynamical processes in 
physics and many other fields of science mathematical 
models are used, like, for example, ordinary differential 
equations (ODEs), partial differential equations (PDEs), 
or iterated maps. Some of these models are derived from 
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first principles while others are the result of a general 
black-box modeling approach (e.g., based on neural net¬ 
works) . These models typically contain two kinds of vari¬ 
ables and parameters: those that can be directly mea¬ 
sured or are known beforehand (e.g., fundamental physi¬ 
cal constants) and others whose values are unknown and 
very difficult to accesd^. To estimate the latter estima¬ 
tion methods have been devised that aim at extracting 
the required information from the dynamics, here rep¬ 
resented by the model equations and the experimentally 
observed dynamical evolution of the underlying process. 
Different approaches for solving this dynamical estima¬ 
tion problem have been devised in the past, including 
(nonlinear) observer or synchronization schemes 3 13 | par¬ 
ticle filter#^, a path integral formalism ! 15 * 16 *, or optimiza¬ 
tion based algorithms 17 "HI. 

Before applying such an estimation method one may 
ask whether the available time series (observable) actu¬ 
ally contains the required information to estimate a par¬ 
ticular unknown value. In control theory this is called 
observability problem and it can for linear systems of 
ODEs be analyzed and a nswe red by means of the so- 
called observability matrix 20 21 . Using derivative coordi¬ 
nates this approac h can be generalized for nonlinear con¬ 
tinuous system^- 0 * 22 * 23 *. For state estimation of chaotic 
systems Letellier, Aguirre and MaquetP®^ considered 
continuous dynamical systems 

x = f (x) (1) 

that generate some observed signal s(t) = h(x(t)) G R 
where x G U C R M is the state of the system, U is a 
smooth submanifold of R M , and h : R M -G R denotes a 
measurement or observation function. Consider now D- 
dimensional derivative coordinated ^ of the observed 
signal s(t) 

y = (s,s,s, ...,s (D_1) ) = F(x) el® (2) 

where s® stands for the k- th temporal derivative of 
s(t), D is the reconstruction dimension , and F is called 
derivative coordinates map 31 . If this map is (at least lo¬ 
cally) invertible, then we can uniquely determine the full 
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state vector x(t) G U from the signal s(t) and its higher 
derivative^! Furthermore, small perturbations in y 
should correspond to small perturbations in x and vice 
versa. Therefore, we want the map F :U F(U) C R D 
to be an immersion , i.e. a smooth map whose deriva¬ 
tive map is one-to-one at every point of U, or equiva¬ 
lently, whose Jacobian matrix DF(x ) has full rank on 
the tangent space (here rank(DF(x)) = M Vx G IX). 
This does not imply that the map F itself is one-to-one 
(.F(x) = F( z) =>* x = z), since the derivative coordi¬ 
nates © may provide states y G R D with two (or more) 
pre-images separated by a finite distance. Therefore, the 
observability analysis presented in the following is local , 
only, because it is based on analyzing (the rank) of the 
Jacobian matrix DF. The (global) one-to-one property 
of the map F is not checked (what would be necessary, 
and for compact U also sufficient, to show that F is an 
embedding ^). 

The D x M Jacobian matrix DF(x) can be computed 
by means of the vector field given in Eq. (flj) . In fact, 
for linear ODEs the Jacobian matrix DF(x) conforms 
with the observability matrix known from (linear) con¬ 
trol theory 2 '-! To es timat e the rank of DF(x) Letellier, 
Aguirre, and Maquet^^U suggest to compute the eigen¬ 
values /i/e > 0 of the M x M - matrix 

A(x) = DF tr (x) • DF(x). (3) 


Nonzero eigenvalues indicate full rank of DF(x) and thus 
local invertibility of F at x. To quantify the (local) in¬ 
vert ibility of F(x) and thus the (local) observability of 
the full state x Aguirre, Letellier, and MaquetP^ 5 intro¬ 
duced the observability index 


(5(x) 


l^min 

l^max 


(A 

(A 


(4) 


where /i m ^ n (A) and /a max (A) denote the smallest and the 
largest eigenvalue of the matrix A , respectively. Time 
averaging (along the available trajectory for 0 < t < T) 
yields 


map) where all model parameters are known and only 
state variables have to be estimated from the observed 
time series. 


A. Estimating state variables of a filtered Henon map 

For an M dimensional discrete system 

x(n + 1) = g (x(n)) (6) 

which generates the times series {s(n)} with s(n) = 
h(x(n)) where n = 1, ..., 7V we can construct D dimen¬ 
sional delay coordinateX9M3\ with reconstructed states 

y (n) = (s(n),s(n + 1),...., s(n + D - 1)) (7) 

= G+(x(n)) G R d . 

Again we assume that all states of interest x he within 
a smooth manifold U C M m . Here we consider delay co¬ 
ordinates forward in time. The function G + is therefore 
called forward delay coordinates map G + : U G+(U ) C 
R d . It is also possible to use delay coordinates backward 
in time, or mixed forward and backward, and we shall 
address this issue in Sec. EE 

As already discussed with derivative coordinates in the 
previous section a state x = (aq,...,%) is locally ob¬ 
servable from the time series {s(n)} if G is an immersion, 
i.e. if the Jacobian matrix DG{x. ) has maximal (full) 
rank M at x. The corresponding Dx M Jacobian matrix 
DG(x) can be computed using the iterated map Eq. 

If the Jacobian matrix DG{x ) has maximal rank M (as¬ 
suming M < D) then G is locally invertible (on G(U)). 
“Local” means that still a delay vector y could possess 
different pre images (separated by a finite distance). 

To motivate and illustrate this analysis we consider the 
Henon map 


xi(n + 1) = 1 — ax\{n) + bx 2 (n) (8) 

x 2 {n + 1) = x\{n) (9) 


S= f J o <5 ( x (^)) di - ( 5 ) 

Instead of derivative coordinates we consider in the 
following delay coordinates 29 33 . Furthermore, we ex¬ 
tend the observability analysis to parameter estimation 
and compute a specific measure of uncertaint}®l for each 
state variable or parameter to be estimated. Last not 
least, we are not only interested in quantifying the av¬ 
erage observability (like 5 in Eq. ©) but also in local 
variations that can be exploited during the state and pa¬ 
rameter estimation process. 

II. DELAY COORDINATES AND OBSERVABILITY 


with parameters a = 1.4 and b = 0.3. In the following we 
shall assume that the dynamics of this system is observed 
via a filtered signal s(n) provided by an FIR-filter 

s(n) = x\(n) + cxi(n — 1) 

= xi(n) + cx 2 (n) = h(x.(ri)) (10) 

with filter parameter c. 

For two dimensional delay coordinates the delay coor¬ 
dinates map reads 

G(x(n)) = (s(ri),s(n + 1)) 

= (xi(n) + cx 2 (n), 

1 — axf(n) + bx 2 {n) + caq(n)) 


To motivate the concepts to be presented in the follow¬ 
ing we shall first consider a discrete time system (iterated 


or 


G(x) = (x\ + cx 2 ,1 — ax\ + bx 2 + cx i). (11) 
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The Jacobian matrix of the map G is given by: 

CG M-(- 2 J 1 + cl) < 12 > 

and its determinant 

det(DG(x)) = 2acxi +b — c 2 (13) 

vanishes for all states x = (xi,x 2 ) on the singular line 


For c —y 0 the FIR-filter is (asymptotically) deactivated 
and the critical line disappears (c 0 => x\ — 00 ). 
For 0.0867 < c < 3.66, however, the critical line crosses 
the chaotic attractor as shown for c = 0.5 in Fig. 0 - 






FIG. 1. (Color online) (a), (b) Henon attractor and singular 
axes T-(short, blue) and (longer, red) for filter 

parameter c = 0.5. (c), (d) Delay coordinates for c — 0.5 
(xf = -0.0357) and c = 0.08 (xf = -1.311). 


B. Forward, backward, and mixed delay coordinates 

Instead of using state space reconstruction based on 
forward delay coordinates 0 one could also use backward 
delay coordinates 

y(n) = (s(n), s(n — 1),, s(n — D + 1)) (15) 

= G_(x(n)) G R d 

or more general, a combination of forward and backward 
components 

y(n) = (s(n - £>_),... s(n — 1), s(n), s(n + 1), 

...,s(n + £>+)) (16) 

= G±(x(n) ;D-,D+) G 


called mixed delay coordinates in the following, with re¬ 
construction dimension D = 1 + D_ + D + . To obtain the 
backward components s(n — k) = ft((x(n — k)) the in¬ 
verse map x(n — 1) = g _1 (x(n)) and its Jacobian matrix 
are required (here we assume that the dynamics is time 
invertible). For discrete time systems (like the Henon 
example) the underlying map © has to be inverted and 
for continuous time systems the inverse of the flow can in 
principle be computed by integrating the system ODEs 
0 backward in time. In both cases, however, problems 
may occur in practice, because an explicit form of the 
inverse map may not exist and backward integration of 
dissipative systems results in diverging solutions and nu¬ 
merical instabilities (for longer integration times). De¬ 
spite these difficulties inclusion of backward components 
turns out to be beneficial for the estimation task as will 
be demonstrated in the following for the Henon examples 
and the Rossler system. 


C. Noisy observations and uncertainty 

At states (xi,X2) with x\ ^ x\ the delay coordinates 
map G is in principle invertible, but the inverse can be 
very susceptible to perturbations in y like measurement 
noise. To quantify the robustness and the sensitivity of 
the inverse with respect to noise we consider the singular 
value decomposition of the Jacobian matrix DG of the 
delay coordinates map 


DG = U - S - V tr (17) 

where S = diag(<7i,..., <tm) is an M x M diagonal matrix 
containing the singular values G\ > 02 > ... > > 0 

and V = (uW,..., u( M )) and V = (v^,..., v^ M )) are 
orthogonal matrices, represented by the column vectors 
u (d £ and E M m , respectively. V tr is the trans¬ 
posed of V coinciding with the inverse E -1 = V tr . Anal¬ 
ogously, U tr = U~ x and the inverse Jacobian matrix 
reads 

DG- 1 = V • S' -1 • U tr (18) 

where S~ x = diag(l/oi,..., 1/<7 m)- Multiplying by V 
from the right we obtain DG~ 1 U = V • S~ x or 

ZXrV™) = —v (m) (m 35s 1,, M). (19) 

This transformation is illustrated in Fig. [2] and it de¬ 
scribes how small perturbations of y in delay reconstruc¬ 
tion space result in deviations from x in the original state 
space. Most relevant for the observability of the (origi¬ 
nal) state x is the length of the longest principal axis 
of the ellipsoid given by the inverse of the smallest sin¬ 
gular value <jm (see Fig. [ 2 ]). Small singular values cor¬ 
respond to directions in state space where it is difficult 
(or even impossible) to locate the true state x given a 
finite precision of the reconstructed state y. For the fil¬ 
tered Henon map we find that the closer the state x is 
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FIG. 2. (Color online) The inverse Jacobian DG~ 1 (y) maps 
perturbations of y in delay reconstruction space to deviations 
from the state x whose magnitudes depend on the direction 
of the perturbation as described by Eq. (19). 


to the critical line given by x\ (14) the more severe is 
this uncertainty. This is illustrated in Fig. [l^,b where at 
some points x the ellipses spanned by the column vec¬ 
tors of the matrix V • S~ x are plotted. Figure [ 3 ] shows 
(color coded) the logarithm of the ratio smallest singular 
value (Jrnin = &m (here: M = 2) divided by the largest 
singular value a max = o\ vs. state variables x\ and x 2 
in a range of coordinates containing the chaotic Henon 
attractor. In Fig. It D = 2 dimensional forward de¬ 
lay 0 is considered where at x\ the smallest singular 
value (Jmin = o'M = 02 vanishes indicating the singu¬ 
larity ( [14] ) illustrated in Fig. [lji,b. If the reconstruction 
dimension D is increased from D = 2 to D = 3 the sin¬ 
gularity disappears as can be seen in Fig. 11 showing the 
ratio (Jmin jo max (color coded) for D = 3 dimensional for¬ 
ward delay coordinates. For comparison, Figs.^d show 
results obtained with mixed delay coordinates (16]) and 
backward delay coordinates ( [l5| ), respectively. The white 
areas in Fig. |3]i correspond to ratios Jmin/omax < 0-01 
indicating poor observability (due to fast divergence of 
backward iterates of the Henon map). Further increase 
of the reconstruction dimension (D = 4 or D = 5) results 
in even larger values of Omin/omax (not shown here). 

To assess the observability on the Henon attractor we 
computed histograms of ratios amin/omax at 10 6 points. 
Figure [4] shows these histograms for the same coordinates 
used to generate the corresponding diagrams in Fig. [3] 
The best results (large ratios) provide mixed delay coor¬ 
dinates (Fig. [ 4 ]:). We speculate that this is due to the 
fact that forward and backward components cover differ¬ 
ent directions in state space (similar to Lyapunov vectors 
corresponding to positive and negative Lyapunov expo¬ 
nents) . 

If the perturbations of the reconstructed state y are 
due to normally distributed measurement noise they can 
be described by a symmetric Gaussian distribution cen¬ 
tered at y 


Q( y) 


exp [-|(y - y ) tr1 (y - y)] 

v /(27r)- D det(I] y ) 


( 20 ) 


where E y = diagfp 2 ,..., p 2 ) = p 2 Id denotes the D x D 
covariance matrix (Id stands for the ZJ-dimensional unit 
matrix) and the standard deviation p quantifies the 
noise amplitude. For (infinitesimally) small perturba¬ 


(a) 


D= 2 


lo s(^) (b) 


D= 3 


lQ g(S;) 




FIG. 3. (Color online) Local observability of the filtered 
Henon map ([8])-(flO]) . Logarithm of the (color coded) ratio 
of the smallest singular value Jmin — cfm (M = 2) divided 
by largest singular value cr macc = a 1 vs. coordinates x\ and 
X 2 for c — 0.5. (a) D = 2 dimensional forward delay co¬ 

ordinates (I 7 L (b) D — 3 dimensional forward delay coordi¬ 
nates 0.(3 D — 3 dimensional mixed delay coordinates ( |16| ) 
(D- — 1 = D+), and (d) D = 3 dimensional backward delay 
coordinates Jl5l. 


(a) D= 2 (b) D= 3 



logKm./ Gnax ) log (J min/Gnax) 


FIG. 4. (Histograms of rations omin/crmax computed at 10 6 
points of the Henon attractor with: (a) D = 2 dimensional 
forward delay coordinates 0, (b) D — 3 dimensional forward 
delay coordinates 0, (c) D = 3 dimensional mixed delay co¬ 
ordinates ( pT| (D- — l — £>-|_) and (d) D — 3 dimensional 
backward delay coordinates (15). Compare corresponding di¬ 
agrams in Fig. [3] 


tions Ay = y — y this distribution is mapped by the 
(pseudo) inverse of the delay coordinates map to the 
(non-symmetrical) distribution 


P(x) 


exp [— | (x — x) tr E a , 1 (x — x)] 

y(27r) M det(E x ) 


( 21 ) 
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centered at x with the inverse covariance matrix 

S” 1 = DG tr • S" 1 • DG = L DG tr • DG (22) 

= \v ■ S 2 ■ V tr . (23) 

P 

The marginal distribution Pj of a given state variable 
Xj centered at Xj is given by 


Pj(Zj) 


1 




(24) 


where the standard deviation pj is given by the square 
root of the diagonal elements of the covariance matrix 


Pj — 


>33 m 


(25) 


Using Eq. (22) the standard deviation of the marginal 


distribution Pj can be written 


Pj = P\J[DG^ ■ DG\ = P J\V • 5-2 • (26) 

and we consider in the following the factor 



= \J[Dtp • DG} 7/ = J\V • 5-2 • (27) 

as our measure of uncertainty when estimating Xj, be¬ 
cause it quantitatively describes how the initial standard 
deviation p is amplified when estimating variable xj 
Figure |5] shows the uncertainties v\ and v 2 vs. (xi 1 x 2 ) 
for different D = 3 dimensional delay coordinates. In 
Fig. [4^i, b results obtained with D = 3 dimensional for¬ 
ward delay coordinates 0 are given. Large uncertainties 
occur mainly in a vertical stripe located near the singular¬ 
ity at x\ (Eq. @) occurring for D = 2. Figures [ 5 J 3 , d,e,f 
show uncertainties of x\ and x 2 obtained with mixed de¬ 
lay coordinates ((c),(d)) and backward delay coordinates 
((e), (f)). For mixed delay coordinates (Fig. [5]qd) ar¬ 
eas with very high uncertainties occur near the origin, 
but along the attractor v\ and take only relatively low 
values. This is also confirmed by the ^-histograms on the 
attractor given in Fig.[6]for the same delay coordinates as 
used in Fig. [4j Again, the mixed delay coordinates turns 
out to be superior to the purely forward or backward co¬ 
ordinates. Furthermore, the dependance of the range of 
uncertainty values on the type of coordinates is different 
for different variables. While the uncertainty v\ of x\ in¬ 
creases when changing from forward to backward delay 
coordinates (Figs. [6^L,e) , the uncertainty of x 2 exhibits 
the opposite trend (Figs. [6]3,f). 


FIG. 5. (Color online) Uncertainty ( |27| of the variables x\ 
and X 2 of the Henon map for c = 0.5 and different D — 3 
dimensional delay coordinates, (a), (b) forward coordinates 
0. (c), (d) mixed coordinates ( |16| ), and (e), (f) backward 
coordinates ( p~5| ). 



log(iq ) log(^ 2 ) 


FIG. 6. Histograms of uncertainties ( |27| of x\ and X 2 on 
the Henon attractor (computed at 10 6 points) for D — 3 and 
c = 0.5 with (a), (b) forward, (c), (d) mixed, and (e), (f) 
backward delay coordinates (comp. Fig. [4]). 


D. State and parameter estimation 

Until now only the state variables x\ and x 2 are con¬ 
sidered as unknowns to be estimated and the parameters 
a and b of the Henon map and c of the FIR filter are 


assumed to be known. We shall now consider the general 
case including unknown parameters p = (pi,... ,pp) G 
M p of the dynamical system and unknown parameters 
q = (</i,..., </q) G of the measurement function 
h(x, q). Let the dynamical model be a M-dimensional 
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discrete 


where: 


or a continuous 


c(n+ 1) = g(x(n), p) 


x = f(x,p) 


dynamical system generating a flow 


it . ttdM 




dM 


(28) 


(29) 


(30) 


with discrete t = n G Z or continuous t G M time. Fur¬ 
thermore, let’s assume that a time series {s(n)} of length 
N is given observed via the measurement function 


( dh(<f> d ~ t ~ (x,p),q) dh(<f) d ~ t ~ (x,p),q) ^ 

dxi ’ ’ ’ Bxm 


A = 


9h(<t> T ~(x,p),q) 
dx i 

<9/i(x, q) 
dxi 

dh((j) T + (x,p),q) 
dxi 


dh(cf) D +' r + (x,p),q) 


dxi 


dh(<j> T ~ (x,p) ,q) 

dx M s 

<9/i(x,q) 

dxM 

a/i(0 T +(x, p),q) 
dx M 


dh(4> D + T + (x,p),q) 


dx M 


s(t) = /&(<£* (x, p),q) 


(31) 


from a trajectory starting at x. 

This provides the D-dimensional delay coordinates 

y = (s(-D-T-), s(-T-), s(0), s(r+),..., s(£>+t+)) 

= G(x,p,q; D_,D + , r_,r+) G (32) 

with D = 1 + D_ + D + . Here the delay coordinates map 
G is considered as a function of: (i) the state x and the 
parameters p of the underlying system, (ii) the parame¬ 
ters q of the measurement function, (iii) the dimension 
parameters D_ and D + , and (iv) the delay times r_ and 
r + in backward and forward direction, respectively. The 
option to use different delay times, r_ and r + for the 
backward and forward iterations is motivated by the fact 
that for dissipative systems backward solutions (j) T ~ (x) 
quickly diverge and therefore a choice r_ < r + may be 
more appropriate. For the same reason D_ has typically 
to be smaller than D + . Since the reconstruction dimen¬ 
sions and the delay times are chosen a priori and are not 
part of the estimation problem they shall not be listed as 
arguments of G to avoid clumsy notation. The Jacobian 
matrix DG(x, p, q) of G has the structure 


B = 


/ dh((f> T — (x,p) ,q) 
1 dp i 


dh(<f> T ~ (x,p),q) 
dpi 

0 

dh((f) T + (x,p),q) 
dpi 


\ dh(<f> D + T + (x,p),q) 
\ dpi 


Cm 


dh(<j> T ~ (x,p),q) 
dqi 

dh(x, q) 

dot 

d/iQ T +(x, p),q) 
dqi 


dh(4> D + T + (x,p),q) 
dqi 


dh(4> d - t - (x,p),q) \ 
dpp » 


dh(<f> T ~ (x,p),q) 
dpp 

0 

dh((f) T + (x,p),q) 
dpp 


a/i(0 jP + T +(x, p),q) / 

dp P / 


/ dh(<f> d - t - (x,p),q) dh{ <f> D ~ T — (x,p) ,q) \ 

1 dqi ' ' ' dq L ' 


dh(<f> T — (x,p) ,q) 
dqL 

<9/i(x,q) 

dqL 

dh((j) T + (x,p),q) 
dq L 


dh(4> D + T + (x,p),q) . 

dq L / 


DG(x,p,q) = (AB,C) 


(33) 


and it can also be written as 
DG = 

( V x h{<t>~ D - T - (x, p), q) • D x 4>~ d - t - (x, p) \7 x h{<t>- 


'(x,p),q )-D p <j> d - t - (x, p) V q h(4> d - t - (x, p), q) \ 


V x h{4>~ T ~ (x, p), q) • D x (j>~ T ~ (x, p) 
Vxft(x,q) 

V x /i(0+(x, p), q) • D x <f>l(x, P) 


V x h(4> T ~ (x, p), q)) • D p (j> T -(x,p) 
0 

V !C /i(0+(x,p),q)) • D p (p T + (x, p) 


V 9 /i( 0“ t - (x, p), q) 
V,li(x, q) 
V 9 /i(0+(x,p),q) 


\ ^xh((f) D+T+ (x, p), q) • D x (f> D+T + (x, p) V*h(0 D + r +(x,p),q) • Z) p 0 D + T +(x,p) V g ft(0 i3 + T +(x, p), q) / 


(34) 


r 


where 


„ 7 / . i dh dh \ , x . . 

v .Mx,q)=( S; ^)(*.1) < 35 > 

v ,Mx,q)=(^,...,^)(x,q ) (36) 


and D x cj. i>*(x, p) and D p ^(x, p) denote the Jacobian ma¬ 
trices of the flow whose elements are derivatives with 
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respect to the state variables x and the parameters p, re¬ 
spectively. For discrete dynamical systems (28) the Jaco- 
bians Z^^x, p) and Dp^x, p) can be computed using 
the chain rule and the recursion schemes 


Dx4> t+1 {x,p) = Ac50‘(x,p),p) ■ D x 4>\yL, p) (37) 

D p (/) t+1 (x, p) = D x g((p t {x, p),p) • D p g{4>\x, p),p) 

+D p( f) t (x, p) (38) 


with D x (f>°(x,p) = I D (D x D unit matrix and 

Dp</>°(x, p) = 0. If backward iterations are required 
(D- > 0 and r_ > 0) similar recursion schemes exist 
based on the inverse map g~ l (providing </> _t ) and its 
Jacobian matrice D x g~ x and D p g~ x . Instead of recur¬ 
sion schemes one may also use symbolic or automatic 
differentiation 3 ^. For continuous systems (29) the re¬ 
quired Jacobians can be obtained by simultaneously solv¬ 
ing l inear ized systems equations as will be discussed in 
Sec. 


IIF 


Inverse maps (D_ >0 and r_ > 0) may be 
computed via backward integration of the ODEs (at least 
for short periods of time before solutions diverge sig¬ 
nificantly). An extension for multivariate time series is 
straightforward. 


E. Parameter estimation for the Henon map 

We shall now extend the discussion to include not only 
state estimation but also parameter estimation. For bet¬ 
ter readability only forward delay coordinates are consid¬ 
ered in the following, but all steps can also be done with 
mixed or backward delay coordinates, of course. We first 
consider the case where b and c are assumed to be known 
and only a has to be estimated. In this case, M = 2 
unknown variables and P = 1 unknown system param¬ 
eter exists (while Q = 0). Therefore, delay coordinates 
with dimension D = 3 or higher will be used. Figure [7] 
shows the ratio of singular values cr min /cr maa , = cf^/cfi vs. 
(aq, X 2 ) in a plane in M 3 given by p\ = a = 1.4 (and fixed 
parameters b = 0.3 and c = 0.5). For reconstruction di¬ 
mension D = 3 the ration cr m ^ n /cr maa , is very small for 
extended subsets (x, p) = (aq,X 2 ,.Pi) of the plane (white 
stripes in Fig. Eh If the delay reconstruction dimension 
is increased to D = 4 (Fig. Eh these regions shrink or 
disappear. If the dimension D is increased furthermore, 
the delay coordinates map is locally invertible in the full 
range of aq and x 2 values shown in Fig. [ 7 ] (results not 
shown here). 

Now we include P 2 = b in the list of quantities to be 
estimated. Figure [8] shows the ratio of singular values 
Cmin/cmax for reconstruction dimensions D = 4 and 
D = 5. For D — 4 curves with very low singular value ra¬ 
tios amin /& max exist crossing the Henon attractor which 
disappear for D = 5. 

Figure [ 9 ] shows the uncertainties zq,..., z/ 4 (Eq. (|27|)) 
for D = 5. As can be seen, the values of uncertainties 
vary strongly in the X 1 -X 2 plane and still some islands 
with rather large uncertainties exist. 



FIG. 7. (Color online) Logarithm of ratio smallest singular 
value (Jmin = cr 3 divided by largest singular value (Jmax — CTi 
vs. x\ and X 2 for the case of M-\-P — 3 unknown (aq, X 2 ,pi = 
a). The diagram shows the plane p± = a = 1.4 in the three di¬ 
mensional estimation space. The other parameters are b — 0.3 
and c = 0.5. Diagrams (a) and (b) show the results ob¬ 
tained with forward delay reconstruction dimensions D — 3 
and D — 4, respectively. 




FIG. 8 . (Color online) Logarithm of the ratio of singular 
values (Jmin/ & max = (Ja/cti vs. X\ and X 2 for the case of 
M + P = 4 unknown quantities (xi,X 2 ,pi = a,P 2 = b). The 
diagrams show the X 1 -X 2 plane at fixed pi = a = 1.4 and 
P 2 — b — 0.3 in the four dimensional estimation space for 
c = 0.5. Diagrams (a) and (b) show the results obtained with 
forward delay reconstruction dimensions D — 4 and D — 5, 
respectively. 


Similar results are obtained if we include the remain¬ 
ing parameters c in the estimation problem. Scan¬ 
ning the two-dimensional X 1 -X 2 subspace (plane) of the 
M + P = 5 dimensional estimation problem for (x, p) = 
(xi,x 2 ,PuP2,q) with fixed p 1 = a = -1.4, p 2 = b = 0.3, 
and q = c = 0.5 indicates (almost) vanishing smallest 
singular values as long as D < 8. With D = 9 dimen¬ 
sional delay coordinates the Jacobian matrix DG(x, p) 
has clearly full rank almost everywhere within the cho¬ 
sen range (aq,^) £ [—1., 1.] x [— 1., 1.] as can be seen in 
Fig.. [TO^i. Figures To]3,c,d,e,f illustrate the uncertainties 
zq,...,z/ 5 (Eq. d27|) of xi,x 2 ,pi,P 2 ,g, respectively. 


F. Continuous dynamical systems 


To compute the Jacobian matrix .DG(x, p) ( |34| ) of the 
delay coordinates map G we have to compute the gradi¬ 
ents (35) and (36) of the observation function s = h(x, q) 
and the Jacobian matrices ^^(x, p) and ^^(x, p) 
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FIG. 9. (Color online) Estimation of uncertainties Vj (j = 
1 = 4)) of variables (xi,X 2 ) and parameters (pi — 

a,p 2 — b) of the Henon map ^ obtained with D — 5 dimen¬ 
sional forward delay coordinates . 





FIG. 10. Estimation of all variables x = (x±,X 2 ) and model 
parameters p = (pi,p 2 ) = (a, 6) of the Henon map H, 
and the measurement function parameter q — c (FIR fil¬ 
ter |To|). The output of the FIR filter is forward embed¬ 
ded in D = 9 dimensions, (a) Ratio of smallest and largest 
singular value, (b)-(f) Uncertainties Vj of state variables 
and parameters in the plane {(xi, £ 2 ,pi,P 2 , q) : ((^ 1 ,^ 2 ) G 
[-!.,!.] x pi = a — 1.4, p 2 = b = 0.3, q = c = 0.5}. 


containing derivatives of the flow (ft generated by the dy¬ 
namical system (29) with respect to variables Xj and pa¬ 
rameters pjj respectively. The M x M-matrix R s ^(x, p) 
can be computed by solving the linearized dynamical 
equations in terms of a matrix ODE 


jY = D x f((/) t (x,p),p) ■ Y (39) 

where <^(x, p) is a solution of Eq. ( [29] ) with initial value 
x and Y is an M x M matrix that is initialized as 
y(0) = Im, where Im denotes the M x M identity ma¬ 
trix. Similarly, the M x P-matrix D p ft(yi, p) is obtained 
as a solution of the matrix odeP 

^z = D x f(x(t), p )■ Z + D p f(x(t), p) (40) 

with Z( 0) = 0. Solvin g ([39] ) and ( |40| ) simultaneously 
with the system ODEs (|29|) we can compute P x (/) r (x), 
P x 0 2r (x), etc. and use these matrices to obtain the Ja¬ 
cobian matrix DG of the delay coordinates map G ( [34] ) . 
For mixed or backward delay coordinates the required 
components can be computed by integrating the system 
ODE and the linearized ODEs backward in time. 


1. The Rossler system 

To demonstrate the observability analysis for continu¬ 
ous systems we follow Aguirre and Letelliei^ and con¬ 
sider the Rossler system 


±1 = -X2 - Xs 

X 2 = xi + ax 2 (41) 

x 3 = b + x 3 (x 1 - c) 

with a = 0.1, b = 0.1, and c = 14. 

Time series of different observables oq, or x 2 , or x 3 are 
considered, all of them consisting of N m 10000 values 
sampled with At = 0.1. Figure [TT| shows the Rossler at¬ 
tractor where color indicates the uncertainty of estimat¬ 
ing the variable x\ (first column), or x 2 (second column), 
or x 3 (third column) using forward delay coordinates. 
The results in the first row are obtained when observ¬ 
ing xi, while the diagrams in the second and third row 
show results for x 2 or x 3 time series. The reconstruction 
dimension equals D = 7 and the delay time is r = 0.5. 
The bright yellow bullet indicates the state with the low¬ 
est uncertainty. This state and the D — 1 = 6 following 
states plotted as thick red bullets underly the time series 
values that are used for the delay reconstruction. They 
span a window in time of length ( D — 1 )r = 6 • 0.5 = 3 
which is about one half of the mean period of the chaotic 
oscillations T « 6. The lowest uncertainties are obtained 
for states whose reconstruction involves trajectory seg¬ 
ments following the vertical ^ 3 -excursion on the attrac¬ 
tor. In contrast, trajectory segments starting from states 
with poor observability (large uncertainty) are located in 
the flat part of the Rossler attractor. Figures and b 
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show that using x\ time series low values of v\ occur on 
parts for the attractor where v<i is high (and vice versa). 
Interestingly, this is not the case for delay reconstructions 


based on x 2 time series as can be seen in Figs. [TT] l and 
e, where low uncertainties of x\ and X 2 occur in similar 
regions on the attractor. 



FIG. 11. (color online) Color coded Rossler attractors where colors of points representing states are given by logarithms of 
uncertainty values v\ in the first column, v<i in the second column, and V 3 in the third column. All results are computed 
using forward delay coordinates. The diagrams in the first row show results obtained based on an reconstruction of a x\ time 
series, the second row using X 2 as an observable, and the third row uncertainties of estimates from X 3 data. The reconstruction 
dimension is D = 7 for all nine diagrams. The bright yellow bullet indicates the state with the lowest v -value (respectively). 
To estimate this state time series values at this state and at D — 1 = 6 subsequent states indicated by dark (red) bullets are 
used for delay coordinates. 


r 


In Fig. [12] distributions of uncertainty values of the 
Rossler system are shown that were obtained along an 
orbit of TV" = 10000 states sampled with At = 0.1. The 
distributions are shown as color coded histograms, esti¬ 
mated from the relative frequency of occurrence of the 
corresponding vj (in %). All diagrams show the depen- 
dance of the histograms on the delay time r chosen for 
forward delay coordinates (horizontal axis). The recon¬ 
struction dimension is for all cases D = 4 and all three 
parameters are assumed to be known (and are not part 
of the estimation task, i.e. P = 0). In the first row 


(Figs. 12r,d,g) estimations are based on a x\ time series 


from the Rossler system, and in rows two and three, x 2 
and £3 time series are used, respectively. The uncertain¬ 
ties Vj of the given observable Xj (Figs. 12r,e,i) mostly 
equal one (log 10 (^) ~ 0 ) or are smaller (due to the ad¬ 
ditional information provided by the delay coordinates). 
In general, lowest uncertainties for all variables are ob¬ 
tained when using x\ time series (Figs. fl2^i,d ,g) while X 3 
data provide highest uncertainties (Figs. | 12 |:,f,i). 


Figure [13] shows the same diagrams but now com¬ 
puted using four dimensional mixed delay coordinates 
with D_ = 1 and D+ = 2 . Similar to the results ob¬ 
tained with the Henon map the uncertainties computed 
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FIG. 12. (Color online) Histograms of uncertainties Vj of the Rossler system (41) vs. delay time r. All parameters are assumed 
to be known (M = 3, P = 0) and forward delay coordinates with dimension D — 4 are used. In the first row a x\ time series 
is given, in the second row X 2 data, and in the third row the delay reconstruction is based on x 3 . The three columns show 
histograms of the (logarithm of the) uncertainties z/i, z^ 2 , and z/ 3 , respectively. 
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FIG. 13. (Color online) Histograms of uncertainties Vj of the Rossler system (41) vs. delay time r. All parameters are assumed 
to be known (M = 3, P — 0) and mixed delay coordinates with dimension D = 1 + D- + D+ •= 4 are used (D- m 1, D+ = 2). 
In the first row axi time series is given, in the second row X 2 data, and in the third row the delay coordinates are based on 
x 3 . The three columns show histograms of the (logarithm of the) uncertainties and z/ 3 , respectively. 





























































































11 


for mixed delay coordinates are typically smaller than 
those obtained with forward coordinates. Furthermore, 
the histograms shown in Fig. [13] suggest that for mixed 
delay coordinates using X 2 as measured variable provides 
the best results, followed by the x\ time series. This is in 
contrast to forward coordinates (Fig. [l2| where x\ data 
yield the smallest uncertainties for the other variables {x 2 
and xs). Similar results have been obtained with three 
dimensional forward or mixed coordinates. 

The fact that mixed delay coordinates provide the low¬ 
est uncertainties when using x 2 time series is consistent 
with results for derivative coordinates obtained by Letel- 
lier et al. 24 who found a ranking x 2 > x\ > X 3 (for a dif¬ 
ferent set of model parameters). For better comparison 
with their results we computed the (attractor) average 



of the ratio 



7(x) = 


<£ <n (£>G(x)) 

<7max( DG (, x )) 


(42) 


(43) 


that provides the delay reconstruction analog 7 of the 
observability index Q. Figure [l4| shows 7 vs. the delay 
time r for different delay coordinates (rows) and different 
measured time series (columns). While for forward delay 
coordinates the largest values of the observability index 
occur if x\ is measured (Fig. [bTfc), X 2 time series provide 
best observability if mixed delay coordinates are used 
(Fig. p|i). Note that in most cases high observability 
occurs for r ~ 1 which is very close to the first zero of the 
autocorrelation function (that is often used as preferred 
value for delay reconstruction). 

Fig. [15] shows similar histograms but now for the full 
estimation problem (M = 3 variables and P = 3 pa¬ 
rameters). Forward delay coordinates are used and the 
reconstruction dimension is increased to D = 13 and an 
x\ time series of length N = 10000 is used (with sam¬ 
pling time At = 0.1). For delay times r that are an 
integer multiple of half of the mean period T/2 « 3 rela¬ 
tively high uncertainties of occur, in particular for ^4, 
z/ 5 , and vq. This is due to the well known fact that for 
these delay times the attractor reconstruction results in 
points scattered near a straight (diagonal) line (an effect 
that also occurs when considering delay coordinates of a 
sinusoidal signal). 


III. CONCLUSION 


FIG. 14. Mean observability indices ( |42| ) of the Rossler system 
( |41| vs. delay time r based on three dimensional (dashed 
lines) and four dimensional (solid lines) delay coordinates. 
Left column ((a), (c), (e) forward delay coordinates. Right 
column ((b), (d), (f)) mixed delay coordinates with D- = 1 
and D+ — 1 (dashed lines) or D+ = 2 (solid lines). The 
measured time series is in the first, second, and third row the 
variable xi, x 2 , and X 3 , respectively. 
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FIG. 15. (color online) Histograms of uncertainties Vj of the 
Rossler system ( [IT] ) vs. delay time r. Both, all state variables 
(M = 3) and all parameters (P = 3) are assumed to be un¬ 
known and have to be estimated from a x\ time series. Results 
are obtained using forward delay coordinates with D — 13. 
The uncertainties ^i,^ 2,^3 correspond to state variables xi, 
X 2 , and X 3 , while 124 , ^ 5 , vq quantify uncertainties of estimated 
parameters pi , P 2 , and P 3 . 


Starting from the question “Does some particular 
(measured) time series provide sufficient information for 
estimating a state variable or a model parameter of inter¬ 
est” we revisited the observability problem for nonlinear 
(chaotic) dynamical systems. In particular we considered 
delay coordinates and the ability to recover not measured 
state variables and parameters from delay vectors. This 


requires to “invert” the delay coordinates construction 
process which is at least locally possible, if the Jaco¬ 
bian matrix of the delay coordinates map has maximum 
(full) rank. Furthermore, we investigated how states near 
the delay vector are mapped back to the state and pa¬ 
rameter space of the systems. In this way it is possible 
to quantify the amplification of small perturbations in 
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delay reconstruction space in different directions of the 
state and parameter space. This reasoning gave rise to 
the concept of uncertainties of estimated variables and 
parameters. Both, observability and uncertainties may 
vary considerably in state space and on a given (chaotic) 
attractor. This feature was demonstrated with a discrete 
time system (filtered Henon map) and a continuous sys¬ 
tem (Rossler system). Local observability and uncertain¬ 
ties also depend on the available measured variable (time 
series) and the type of delay coordinates. Best results 
were obtained with mixed delay coordinates, containing 
at least a one step backward in time. 

The obtained information about (local) uncertainties 
in state and parameter estimation can be used in several 
ways for subsequent analysis. First of all, it may help 
to decide whether the planned estimation task is feasible 
at all or whether another observable has to be measured 
instead. For continuous time systems relevant time scales 
(delay times) can be identified where uncertainties are 
minimal. 

The strong variations of local uncertainty values in 
state space (along a trajectory) occurring with the ex¬ 
amples shown here are typical and should be taken into 
account by any estimation method. If the system is in a 
state where, for example, the uncertainty v\ of the first 
variable is high then it might be better not to try to esti¬ 
mate this variable in this state or close to it, because the 
estimate might be poor and may spoil the overall results. 
Instead it makes more sense to wait until the trajectory 
enters a region of state space where x\ can be estimated 
more reliably from the given time series. 

The concrete implementation of such an adaptive ap¬ 
proach depends on details of the estimation algorithm. 
For Newton-like algorithms, for example, it may consist 
of a simple strategy decreasing correction step sizes. An¬ 
other potential application of uncertainty analysis is the 
identification of redundant parameters, i.e., parameter 
combinations that provide the same dynamical output. 
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